#include "gabor.h"
#include "gabor_kernel.cu"
#include "cutil_inline.h"

Gabor::Gabor(int width, int height,int nLevels)
{
	//Image width should be aligned to maximum coalesced read/write size
	//for best global memory performance in both row and column filter.
	//Reduce problem size to have reasonable emulation time
	DATA_W_original = width;
	DATA_W_original_16multi = iAlignUpSub(DATA_W_original, 16);
	DATA_H_original = height;
	DATA_W = iAlignUpSub(DATA_W_original_16multi + 2 * KERNEL_RADIUS_PAD, 16);
	DATA_H = DATA_H_original + 2 * KERNEL_RADIUS_PAD;            //size after padding
	KERNEL_RADIUS_PADto16 = DATA_W_original_16multi - DATA_W_original;
	DATA_W_SHIFT = iAlignUpSub(DATA_W + DATA_H - 1, 16);    //w_size after shift once, multiple of 16, for second column filter
	DATA_W_SHIFT2 = iAlignUpSub(DATA_W_SHIFT + DATA_H - 1, 16);   //w_size after shift once, multiple of 16, for combine and demodulation

	DATA_SIZE_shift1 = DATA_W_SHIFT * DATA_H * sizeof(float);
	DATA_SIZE_shift2 = DATA_W_SHIFT2 * DATA_H * sizeof(float);

	DATA_SIZE_ori_f = DATA_W_original * DATA_H_original * sizeof(float);
	DATA_SIZE_ori_c = DATA_W_original * DATA_H_original * sizeof(Complexgf);
	DATA_SIZE_pad_f = DATA_W * DATA_H * sizeof(float);
	DATA_SIZE_pad_c = DATA_W * DATA_H * sizeof(Complexgf);


	//parameter for necessary block computing
	//calculate parameter for necessary block computing

	filter_interval = iDivUpSub(DATA_W + COLUMN_TILE_H - 1, COLUMN_TILE_W);
	shift_interval = iDivUpSub(DATA_W + BLOCK_DIM - 1, BLOCK_DIM);
	filter_interval2 = iDivUpSub(DATA_W + COLUMN_TILE_H - 1, COLUMN_TILE_W) + 1;
	shift_interval2 = iDivUpSub(DATA_W + BLOCK_DIM - 1, BLOCK_DIM) + 1;
	filter_cof = iDivUpSub(COLUMN_TILE_H, COLUMN_TILE_W);

	
	//set up filter execution parameter
	
	blockGridRows.x = iDivUpSub(DATA_W, ROW_TILE_W);
	blockGridRows.y = DATA_H;

	blockGridColumns.x = iDivUpSub(DATA_W, COLUMN_TILE_W);
	blockGridColumns.y = iDivUpSub(DATA_H, COLUMN_TILE_H);

	blockGridColumns2.x = iDivUpSub(DATA_W_SHIFT, COLUMN_TILE_W);
	blockGridColumns2.y = iDivUpSub(DATA_H, COLUMN_TILE_H);

	threadBlockRows.x = KERNEL_RADIUS_ALIGNED_X + ROW_TILE_W + KERNEL_RADIUS_X;

	threadBlockColumns.x = COLUMN_TILE_W;
	threadBlockColumns.y = 8;

	
	//set up padding execution paramter
	grid_ori.x = iDivUpSub(DATA_W_original_16multi, BLOCK_DIM);
	grid_ori.y = iDivUpSub(DATA_H_original, BLOCK_DIM);
	
	//for the image levels
	grid_oriND.x = grid_ori.x;
	grid_oriND.y = grid_ori.y;
	grid_oriND.z=nLevels;

	// setup shift, modulation, demodulation execution parameters

	grid.x = iDivUpSub(DATA_W, BLOCK_DIM);
	grid.y = iDivUpSub(DATA_H, BLOCK_DIM);
	
	// for the image levels
	gridND.x =  grid.x;
	gridND.y =  grid.y; 
	gridND.z=nLevels;
	
	grid_shift1.x = iDivUpSub(DATA_W_SHIFT, BLOCK_DIM);
	grid_shift1.y = iDivUpSub(DATA_H, BLOCK_DIM);
	
	grid_shift2.x = iDivUpSub(DATA_W_SHIFT2, BLOCK_DIM);
	grid_shift2.y = iDivUpSub(DATA_H, BLOCK_DIM);
	
	threads.x = BLOCK_DIM;
	threads.y = BLOCK_DIM;
	

	//allocate memory

    // memory allocation for gabor computation    

	cutilSafeCall( cudaMalloc( (void **)&d_DataA_c, DATA_W*DATA_H*sizeof(float)) );
	cutilSafeCall( cudaMalloc( (void **)&d_DataA_s, DATA_W*DATA_H*sizeof(float)) );
	cutilSafeCall( cudaMalloc( (void **)&d_Temp, DATA_W*DATA_H*sizeof(float)) );

    cutilSafeCall( cudaMalloc( (void **)&d_DataA_padto16multi, DATA_W_original_16multi*DATA_H_original*sizeof(float)) );
	cutilSafeCall( cudaMalloc( (void **)&d_DataA_real, DATA_W*DATA_H*sizeof(float)) );
	cutilSafeCall( cudaMalloc( (void **)&d_DataA_imag, DATA_W*DATA_H*sizeof(float)) );
    cutilSafeCall( cudaMalloc( (void **)&d_DataB, DATA_W*DATA_H*sizeof(float)) );
    cutilSafeCall( cudaMalloc( (void **)&d_shift1, DATA_SIZE_shift1) );
    cutilSafeCall( cudaMalloc( (void **)&d_Temp_shift, DATA_SIZE_shift1) );
    cutilSafeCall( cudaMalloc( (void **)&d_shift2, DATA_SIZE_shift2) );

}

Gabor::~Gabor()
{
	cutilSafeCall( cudaFree(d_DataA_c) );
	cutilSafeCall( cudaFree(d_DataA_s) );
	cutilSafeCall( cudaFree(d_Temp) );
	cutilSafeCall( cudaFree(d_DataA_padto16multi) );
	cutilSafeCall( cudaFree(d_DataA_real) );
	cutilSafeCall( cudaFree(d_DataA_imag) );
	cutilSafeCall( cudaFree(d_DataB) );
	cutilSafeCall( cudaFree(d_shift1) );
	cutilSafeCall( cudaFree(d_Temp_shift) );
	cutilSafeCall( cudaFree(d_shift2) );
	
}


//odata, idata are both device memory, C language interface
void Gabor::pad_GPU(float* odata, float* idata)
{
	if(KERNEL_RADIUS_PADto16 == 0)
	{
		CUDA_SAFE_CALL( cudaMemcpy(d_DataA_padto16multi, idata, DATA_W_original_16multi * DATA_H_original * sizeof(float), cudaMemcpyDeviceToDevice) ); 
	}
	else
	{
		padtomultiple16<<<grid_ori, threads>>>(
		d_DataA_padto16multi,
		idata,
		DATA_W_original_16multi,
		DATA_H_original,
		KERNEL_RADIUS_PADto16
		);
	}

	//pad image
	padinside<<<grid_ori, threads>>>(
		odata,
		d_DataA_padto16multi,
		DATA_W_original_16multi,
		DATA_H_original,
		KERNEL_RADIUS_PAD
		);
		
	padleftright<<<grid, threads>>>(
		odata,
		DATA_W,
		DATA_H,
		KERNEL_RADIUS_PAD
		);
		
	padupdown<<<grid, threads>>>(
		odata,
		DATA_W,
		DATA_H,
		KERNEL_RADIUS_PAD
		);
}
void Gabor::pad_GPU_ND(float* odata, float* idata)
{
    if(KERNEL_RADIUS_PADto16 == 0)
	{
		CUDA_SAFE_CALL( cudaMemcpy(d_DataA_padto16multi, idata, DATA_W_original_16multi * DATA_H_original * sizeof(float), cudaMemcpyDeviceToDevice) ); 
	}
	else
	{
	    padtomultiple16ND<<<grid_oriND, threads>>>(
		d_DataA_padto16multi,
		idata,
		DATA_W_original_16multi,
		DATA_H_original,
		KERNEL_RADIUS_PADto16
		);
	}
	
	padinsideND<<<grid_oriND, threads>>>(
	odata,
	d_DataA_padto16multi,
	DATA_W_original_16multi,
	DATA_H_original,
	KERNEL_RADIUS_PAD
	);
		
	padleftrightND<<<gridND, threads>>>(
	odata,
	DATA_W,
	DATA_H,
	KERNEL_RADIUS_PAD
	);
	
	padupdownND<<<gridND, threads>>>(
		odata,
		DATA_W,
		DATA_H,
		KERNEL_RADIUS_PAD
		);

}

void Gabor::center_surround_GPU(float* odata, float* idata, Filter_kernel& kernel_c, Filter_kernel& kernel_s)
{
	gaussian_GPU(d_DataA_c, idata, kernel_c);
	gaussian_GPU(d_DataA_s, idata, kernel_s);

	minus<<<grid, threads>>>(
		odata,
		d_DataA_c,
		d_DataA_s,
		DATA_W,
		DATA_H
		);
}
void Gabor::center_surround_GPU_ND(float* odata, float* idata, Filter_kernel& kernel_c, Filter_kernel& kernel_s)
{
	gaussian_GPU(d_DataA_c, idata, kernel_c);
	gaussian_GPU(d_DataA_s, idata, kernel_s);

	minusND<<<gridND, threads>>>(
		odata,
		d_DataA_c,
		d_DataA_s,
		DATA_W,
		DATA_H
		);
}

//idata, odata are both device memory, h_kernel_x, h_kernel_y, h_kernel_t are all host kernel, x, y, diagonal axis respectively, C language interface
void Gabor::gabor_GPU(Complexgf* odata, float* idata, Filter_kernel& kernel)
{
	//modulation
	modulation<<<grid, threads>>>(
			d_DataA_real,
			d_DataA_imag,
		    idata,
			kernel.omegax,
			kernel.omegay,
			DATA_W,
			DATA_H
		);
	CUT_CHECK_ERROR("modulation execution failed\n");
	
	//anigaussian
	gaussian_GPU(d_DataA_real, d_DataA_real, kernel);
	gaussian_GPU(d_DataA_imag, d_DataA_imag, kernel);

	//combine and demodulation
	combine<<<grid, threads>>>(odata, d_DataA_real, d_DataA_imag, DATA_W, DATA_H);
	CUT_CHECK_ERROR("combine execution failed\n");

	demodulation<<<grid, threads>>>(odata, kernel.omegax, kernel.omegay, DATA_W, DATA_H);
	CUT_CHECK_ERROR("demodulation execution failed\n");
}


// cut function because the image is padded first for the consideration of border condition, C language interface
void Gabor::cut_GPU(Complexgf* odata, Complexgf* idata)
{
	cut<<<grid, threads>>>(odata, idata, DATA_W, DATA_H, DATA_W_original, DATA_H_original, KERNEL_RADIUS_PAD, KERNEL_RADIUS_PAD);
}

void Gabor::andImage_GPU(Complexgf* odata, Complexgf* idata,Complexgf* imask)
{
    //incomplete
    //andImage<<<grid, threads>>>(odata, idata,imask,DATA_W_original);
}
//Gaussian function
//odata: output
//idata: input
//h_kernel_x, h_kernel_y, h_kernel_t: kernel in the x, y, t axis
//tanp: 0, +1, -1, means the orientation of the gaussian fitler(the angle between the major axis of gaussian filter and x axis)
//      0: phi = 0, or phi = pi/2
//		+1: 0<phi<pi/2
//      -1: pi/2<phi<pi
//for the circular gaussian, you just need the kernel of x and y axis. So the default value of h_kernel_t and tanp is NULL and 0 respectively.
//for the anisotropic gaussian, you need use the cal_parameter_gaussian function to caculate the paramters necessary.
void Gabor::gaussian_GPU(float* odata, float* idata, Filter_kernel& kernel)
{
	cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_X, kernel.h_kernel_x, KERNEL_SIZE_X) );
    cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_Y, kernel.h_kernel_y, KERNEL_SIZE_Y) );
   

	convolutionRowGPU<<<blockGridRows, threadBlockRows>>>(
            d_DataB,
            idata,
            DATA_W,
            DATA_H
    );
    CUT_CHECK_ERROR("convolutionRowGPU() execution failed\n");

    convolutionColumnGPU<<<blockGridColumns, threadBlockColumns>>>(
           odata,
           d_DataB,
           DATA_W,
           DATA_H,
           COLUMN_TILE_W * threadBlockColumns.y,
           DATA_W * threadBlockColumns.y
    );
	CUT_CHECK_ERROR("convolutionColumnGPU() execution failed\n");

	if(kernel.tanp == 0)
	{
		;
	}
	else if(kernel.tanp > 0)
	{
		cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_T, kernel.h_kernel_t, KERNEL_SIZE_T) );
	   shift45<<< grid, threads >>>(d_shift1, odata, DATA_W, DATA_H, DATA_W_SHIFT);   
       CUT_CHECK_ERROR("shift45() execution failed\n");
            
       //columnfilter == diagonal filter
       convolutionColumnGPU45<<<blockGridColumns2, threadBlockColumns>>>(
            d_Temp_shift,
            d_shift1,
            DATA_W_SHIFT,
            DATA_H,
            COLUMN_TILE_W * threadBlockColumns.y,
            DATA_W_SHIFT * threadBlockColumns.y,
            filter_interval,
            filter_cof
        );
        CUT_CHECK_ERROR("convolutionColumnGPU() execution failed\n");
        
        //shift back
        shift135choosing<<< grid_shift1, threads >>>(d_shift2, d_Temp_shift, DATA_W_SHIFT, DATA_H, DATA_W_SHIFT2, shift_interval);   
        CUT_CHECK_ERROR("shift135choosing() execution failed\n");
		
		//cut to make sure the size of input is the same as input
		cut<<< grid, threads >>>(odata, d_shift2, DATA_W_SHIFT2, DATA_H, DATA_W, DATA_H, DATA_H - 1, 0);
		CUT_CHECK_ERROR("cutonece execution failed\n");
	}
	else
	{
		 cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_T, kernel.h_kernel_t, KERNEL_SIZE_T) ); 
		// first shift for 135 degree diagonal filter
		shift135<<< grid, threads >>>(d_shift1, odata, DATA_W, DATA_H, DATA_W_SHIFT);   
        CUT_CHECK_ERROR("shift135() execution failed\n");
            
        //columnfilter == diagonal filter
        convolutionColumnGPU135<<<blockGridColumns2, threadBlockColumns>>>(
            d_Temp_shift,
            d_shift1,
            DATA_W_SHIFT,
            DATA_H,
            COLUMN_TILE_W * threadBlockColumns.y,
            DATA_W_SHIFT * threadBlockColumns.y,
            filter_interval2,
            filter_cof
        );
        CUT_CHECK_ERROR("convolutionColumnGPU() execution failed\n");
        
        //shift back
        shift45choosing<<< grid_shift1, threads >>>(d_shift2, d_Temp_shift, DATA_W_SHIFT, DATA_H, DATA_W_SHIFT2, shift_interval2);   
        CUT_CHECK_ERROR("shift45choosing() execution failed\n");

		//cut to make sure the size of input is the same as input
		cut<<< grid, threads >>>(odata, d_shift2, DATA_W_SHIFT2, DATA_H, DATA_W, DATA_H, DATA_H - 1, 0);
		CUT_CHECK_ERROR("cutonece execution failed\n");
	}
}

void Gabor::gaussian_GPU_ND(float* odata, float* idata, Filter_kernel& kernel)
{
	cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_X, kernel.h_kernel_x, KERNEL_SIZE_X) );
    cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_Y, kernel.h_kernel_y, KERNEL_SIZE_Y) );
   

	convolutionRowGPU<<<blockGridRows, threadBlockRows>>>(
            d_DataB,
            idata,
            DATA_W,
            DATA_H
    );
    CUT_CHECK_ERROR("convolutionRowGPU() execution failed\n");

    convolutionColumnGPU<<<blockGridColumns, threadBlockColumns>>>(
           odata,
           d_DataB,
           DATA_W,
           DATA_H,
           COLUMN_TILE_W * threadBlockColumns.y,
           DATA_W * threadBlockColumns.y
    );
	CUT_CHECK_ERROR("convolutionColumnGPU() execution failed\n");

	if(kernel.tanp == 0)
	{
		;
	}
	else if(kernel.tanp > 0)
	{
		cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_T, kernel.h_kernel_t, KERNEL_SIZE_T) );
	   shift45<<< grid, threads >>>(d_shift1, odata, DATA_W, DATA_H, DATA_W_SHIFT);   
       CUT_CHECK_ERROR("shift45() execution failed\n");
            
       //columnfilter == diagonal filter
       convolutionColumnGPU45<<<blockGridColumns2, threadBlockColumns>>>(
            d_Temp_shift,
            d_shift1,
            DATA_W_SHIFT,
            DATA_H,
            COLUMN_TILE_W * threadBlockColumns.y,
            DATA_W_SHIFT * threadBlockColumns.y,
            filter_interval,
            filter_cof
        );
        CUT_CHECK_ERROR("convolutionColumnGPU() execution failed\n");
        
        //shift back
        shift135choosing<<< grid_shift1, threads >>>(d_shift2, d_Temp_shift, DATA_W_SHIFT, DATA_H, DATA_W_SHIFT2, shift_interval);   
        CUT_CHECK_ERROR("shift135choosing() execution failed\n");
		
		//cut to make sure the size of input is the same as input
		cut<<< grid, threads >>>(odata, d_shift2, DATA_W_SHIFT2, DATA_H, DATA_W, DATA_H, DATA_H - 1, 0);
		CUT_CHECK_ERROR("cutonece execution failed\n");
	}
	else
	{
		 cutilSafeCall( cudaMemcpyToSymbol(d_Kernel_T, kernel.h_kernel_t, KERNEL_SIZE_T) ); 
		// first shift for 135 degree diagonal filter
		shift135<<< grid, threads >>>(d_shift1, odata, DATA_W, DATA_H, DATA_W_SHIFT);   
        CUT_CHECK_ERROR("shift135() execution failed\n");
            
        //columnfilter == diagonal filter
        convolutionColumnGPU135<<<blockGridColumns2, threadBlockColumns>>>(
            d_Temp_shift,
            d_shift1,
            DATA_W_SHIFT,
            DATA_H,
            COLUMN_TILE_W * threadBlockColumns.y,
            DATA_W_SHIFT * threadBlockColumns.y,
            filter_interval2,
            filter_cof
        );
        CUT_CHECK_ERROR("convolutionColumnGPU() execution failed\n");
        
        //shift back
        shift45choosing<<< grid_shift1, threads >>>(d_shift2, d_Temp_shift, DATA_W_SHIFT, DATA_H, DATA_W_SHIFT2, shift_interval2);   
        CUT_CHECK_ERROR("shift45choosing() execution failed\n");

		//cut to make sure the size of input is the same as input
		cut<<< grid, threads >>>(odata, d_shift2, DATA_W_SHIFT2, DATA_H, DATA_W, DATA_H, DATA_H - 1, 0);
		CUT_CHECK_ERROR("cutonece execution failed\n");
	}
}


